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1. introduction 

As the molecule that carries all genetic information, DNA has been of great interest 
to chemists and physicists since before the unlocking of its structure by Watson and 
Crick pQ, and both the structure and dynamics of the molecule have been investigated 
in great detail. Another key feature of DNA is the ease of its assembly as a linear 
macro-molecule of truly macroscopic scale; a potential realisation of a quantum wire (if 
it were a conductor). One can also use the DNA molecule as a scaffold to grow quantum 
wires, e.g., by deposition of conducting molecules in the minor groove. [21 El HI 13 

With the experiments by Dekker and his group [3 El Ej on the conduction of single 
strands of DNA, and the controversy about the interpretation of the experiments, which 
suggested that DNA might be a conductor, a veritable floodgate for theoretical studies 
has opened. Most of these consider the conduction (or charge transport) along a one- 
dimensional "pi-stack" (a channel built from overlapping pi orbitals) in the centre of 
the DNA molecule. 

Disregarding electronic transport, mathematical modelling of DNA has a long 
tradition, see for example the book Ref. [3- In recent years the Peyrard-Bishop 
model (a simple two-dimensional model for DNA, which has been used successfully 
to describe the denaturing transition of that molecule) has become quite popular, see 
the review in ^3 for a more complete set of references. Mingaleev et al ^21 Ej have 
produced models which allow them to perform mathematical studies of DNA and its 
properties. This is an active field, and the number of relevant publications is large. The 
most relevant literature can be traced from a number of recent references, e.g., Refs. 
[13 El 113 Hi El CHI EH 123 12U 122 123 and the review [23- 

Recently one of the present authors (together with Dr. Hartmann) [23 looked at the 
model of Mingaleev et al and proposed a further generalisation of their one chain model 
to two chains. This two chain generalisation involved similar terms in the Lagrangian 
(for each chain) to those of the model of Mingaleev et al plus an extra term involving 
the interaction between the chains. The obtained results were encouraging, although 
due to the specific choice of terms, the numerical simulations were very slow. 

In this note we reinvestigate these models and propose a new generalisation which, 
at least, has the advantage of fast numerical convergence. We use this model to 
investigate various aspects of DNA-like substances. Our model can be compared to the 
model used by Bishop et al. [22] , but differs in many aspects. The two most crucial ones 
are the inclusion of the helicity of the molecule and the fact that in our model electrons 
are (correctly or incorrectly, that remains to be seen) allowed to propagate along the 
backbone rather than along the "pi-stack" of the overlapping molecular orbitals of the 
base pairs. The latter assumption is already quite old, and was first proposed by Eley 
[23 m 1962. In this work we concentrate our attention on the mathematical aspects 
of the problem postponing detailed aspects of the applicability of the model to real, 
physical and biological, systems to a further publication. 

In the next section we say a few words about the DNA systems and in the following 
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Figure 1. A schematic model of B-DNA 

section we describe our ideas of how to build a model for DNA-inspired structures which 
can be used to describe such systems. 

In the final section we describe in detail our first studies involving this model. 

2. Some aspects of DNA 

The basic structure of DNA varies quite considerably depending on its environment. B- 
DNA, the standard double helix, actually only occurs in aqueous solution. An idealised 
structure of DNA is given in figure 

We can easily see that each base pair advances the rotation of the chain by about 
36°. The difference between the major and minor grooves is caused by the fact that the 
base-pairs are not exactly at right angles to the backbone. There is also a "propeller 
twist": the plane of the base pairs is twisted like a propeller. See Ref. [27] for more 
details. As can be seen in the figure, the different base pairs are quite close in the centre 
of the molecule. This is why it is believed that DNA gains rigidity by overlapping pi- 
orbitals between the different base pairs. This same mechanism has been invoked in 
allowing charge transport through DNA. 
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Figure 2. A simple model of a section of straight DNA 
3. Modelling DNA 

Our DNA model will consist of two parts: a simple classical model of a chiral molecule, 
that can twist and be locally compressed, but where we disregard the bending modes. 
This would be a natural approach to a stretched DNA molecule, attached between two 
electrical contacts. The situation is sketched schematically in figure El1f- In that figure 
the thin (red) lines denote denote the hydrogen bonds, and are modelled to be at straight 
angles to the "virtual backbone" denoted by the dashed line. The thick (blue and green) 
lines denote the two DNA chains. 

At the same time we assume electrons hop from a base pair to a base pair, 
occasionally moving from one side to the other. Our model for this process will be 
based on the discrete non-linear Schrodinger equation, which can be seen as the semi- 
classical limit of a tight-binding model. 

We shall investigate the electronic model first, before concentrating on the molecular 
dynamics. 

3.1. Electron dynamics 

We assume a simple discrete non-linear Schrodinger equation for the movement of the 
electrons through DNA. We model the electron as hopping along sites labelled i, the 
position of a base on the back-bone. We have two electron fields, one on each ribbon, 
which we call ip and (connected to the two sets of coordinates r and s described 

% This may not be the best approach; stretching of the "rungs" may well be ignored at low temperatures, 
but the bending of the backbone seems to be a low-energy degree of freedom. For simplicity we do not 
model the very complicated bending and kinking |28j . 
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above) . 

The action for the electrons, without coupling between the two chains, is taken to 

be 



A= J dt I -hJ2^i,t ~ h-J2& [ ^ t+ 
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This can, for instance, be derived as the semi-classical limit of a tight-binding model, 
with a density- squared term added. Such non-linearities are crucial, since they limit 
the amount of charge that can collect at any one site. It is not obvious that such non- 
linearities are restricted to the on-site interaction, and we could also study the question 
"how different would the dynamics of the system be had we added an electronic term of 
the density-density correlation form J2ij^iXij l^if \4>j\ 2 7" We have no answer to this 
question - and we leave its resolution to some future work. In this paper we focus our 
attention on the study of the system essentially described by (0) together with a further 
potential describing the interchain interaction. As we shall see this system is already 
sufficiently complicated to exhibit many interesting properties. 



3.2. The interaction between chains 

Of course, one of the key elements of the problem is the interaction between the chains 
(and we limit ourselves to consider only the interaction between the two nearest sites 
on the backbones, i.e. i = j). The most naive choice of such an interaction is to take 

v^ = Y,K{^i+^m ■ (2) 

i 

This form of the interchain interaction term is actually rather uninteresting since it leads 
to a simple oscillation of charges from one strand to the other, while their combination 
behaves as a single unit, and we would thus reproduce the usual results of the discrete 
nonlinear Schrodinger Equation (DNLS) + . 

This is due to the fact that Eq. (J2J) does not take into account the chirality of 
the molecules. A possible improvement would involve making the coupling orientation 
dependent: 

V int = K(^m + e- i2 ^*0,) , (3) 

i 

(where g is yet another parameter). Here 9i is the angle of the i sites on the backbones 
relative to some arbitrary direction (see later). This is not satisfactory, however, since 

+ For a good description of properties of this model see e.g. |29| . 
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for a flat molecule the tunnelling would now depend on its orientation! We can think of 
two further improvements: 

First of all we can subtract the average orientation, 

V int = K(e i2 ^-^m + e-^-^m , (4) 

i 

which leads to a rather non-local coupling. This rather artificial approach will not be 
further investigated in this work, since it is not related to the chirality of the DNA 
molecule. 

Secondly, we can invert the transformation of fields obtained from going from Eq. (J2J) 
to Eq. (j3J), which corresponds to a gauge-like transformation of the naive interaction, 
Eq. 0. This implies that we should make the replacement 

V< -> e^, <f>i - ^' ig6i <t>i (5) 
everywhere (i.e., also in all the interaction terms within each one of the backbones). 
For a flat molecule this has no consequence, and the orientation dependence of Eq. 
disappears, but for a curved molecule, as we shall show below, this indeed includes 
effects due to chirality. 

Our modification then implies that we can use Eq. (j2J), but must replace the non- 
local interaction along the backbone by 

£ Jir-i (VWi + Mi) ^ E J *-i P"*^ + e-Wt-Wfflj) .(6) 
This is obviously rotationally invariant. 



3.3. Continuum limit and chirality 
Let us look at the continuum limit of 

Y^J^j^-^j (7) 

and 

Y^-l, ur " "-o]o r (8) 

in order to see that this expression really describes a chiral molecule. The second term is 
simply related to the first by the replacement d{ — > — 0, and ip — > <f>. We thus concentrate 
on the first term, where we assume that Jui falls sufficiently fast so that we need consider 
only the nearest-neighbour terms. Calling the contribution of these terms, divided by 
Ji, k\, we find that 

h = J2 (^i+ie i9(ei+1 ~ dl) + cc.) . (9) 
% 

Next we assume that 9i + \ — Qi is small and we expand k\ to first order in this quantity. 
Before doing so, we introduce some time-saving notation, 

m = W^-ia ^ {fh = W + /,-i /2 - (10) 
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fc i = (^+i(l + ig(69) i+1/2 ) + c.c.) 

i 

= \Y1 + ipi-i + c - c -) 

i 

+ | (^i+l(^)i+l/2 + ^-l(^)i-l/2) + CC.) 

2 



= tEW(^)* +c - c -) + Ei^ 



+ E (^) < + 1 /2(^)<+i/2 - (ii) 1-1/2(88)1-1/2) + c.c.) 

i 

+^ 2 E(^(w)(i|^)]) i +c.c). (ii) 

i 

If we combine the "onsite" term |-?/>j| 2 in a similar term already present in the model, we 
find that the remaining terms become 

h/h =\J dx (rd x (d x + igA)ij + ibd x {d x - igA)^) 

» - J dx {{D X ^)*D X ^) . (12) 

Here A = d x 9, D x = d x + igA In the last step we have added a second order term in 9 
(neglected in our initial evaluation). Thus we end up with a covariant derivative. The 
other chain (0) contains the complex conjugate covariant derivative, 

kt~-Jdx {{D x 4>yD x <t>) . (13) 

Thus we see that with this structure of the coupling to a 'magnetic flux' caused by the 
warping of the molecule through each plaquette, we have indeed constructed a chiral 
model of electrons moving through a molecule. 



3.4- Solitons in a curved background 

First we shall study the discrete solitons (breathers) of the model, i.e., solutions of 
the form ^ = e in '/-M> 4>i — e in */i,i- We assume that the distance between the 
molecules is specified by the parameters d = 2.39 nm, Sz = 3.38/10 = 0.338 nm, 
59 = 2tt/10 = 0.628. 

This leads to the non-linear eigenvalue problem 

2tuvf kti -» :;,A " " 2 -/, .„/a., - X\h,i\ 2 h,i + Kf u = Klfc , (14) 
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amplitude 

Figure 3. The profile of the DNLS breather for various values of the coupling constant 
Jo, starting from Jo = 0.22 and a solution localized at i = 25, increasing Jo in syteps 
of 0.02. The transition occurs between J = 0.3 and Jo = 0.32. Wider curves represent 
larger values of Jo . 
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Figure 4. The transition between solitonic and delocalised solutions as a function of 
the coupling constant Jo- 

where the term containing A = —A denotes the coupling to the other chain. One of the 
ways this problem can be solved is through a constrained minimization, with Ml being 
the Lagrange multiplier imposing charge normalisation. 

We have studied the nature of such solutions, for fko = 1, Ji-j = Joe - ^- 7 , x = 0.1 
and k = 1 We have found that the coupling between the two chains leads to solitons 
that are distributed evenly over the two backbones, but where the phase varies. 

If we ignore the coupling between the two strands of the double helix for a while, 
and analyse the motion of a single electron along an equilibrium chain, we are essentially 
solving the DNLS, and the interesting solutions are breathers of the form <pi(t) = e luJt fi 
• These can be found quite easily for the current model by a simple steepest descent 
method (minimizing the electronic energy). In figure El we present the profiles of the 
DNLS breather for various values of the coupling constant Jo- in our simulations we 
have found a sharp transition from a localised to delocalised solution for values of Jo 
about 0.0282. As can be seen from figure 0] the transition appears to be of first order, 
and the solutions are metastable beyond the transition. 

When we increase the coupling we find that it is instructive to look in more detail 
at the energetics of the problem, as shown in figure where (3 = lnm -1 and various 
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Figure 5. The energy (mass) of the breather for two coupled chains, as a function of 
the coupling constant Jo, and {3 = lnm - . From bottom to top, the curves correspond 
to: black g — 0, red g = 0.5, green g = 1, blue g = 1.2 and violet g = 1.5. The 
dashed (orange) curves show the single chain solution from above. In each case we 
have performed our calculations for both the increasing and decreasing Jo. 
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Figure 6. The profile of the DNLS breathers as a function of the coupling constant 
of the coupling constant Jo, and (3 = lnm -1 . The left figure shows g = 0, the right 
one g = 1.5. The sharpest curves occur for the smallest value of Jo, with a sequential 
broadening on increase. 

values of g have been used. We still see the phase trasition, but note that there is a shift 
in the position, and that it seems to change from first to second or higher order with 
increasing chiral charge g. As can be seen in figure EH there is a change in the nature of 
the solitonic positions as well, and we are no longer able to find fully delocalised ones. 
If we increase j3 the transition seems to weaken even more, and the solitons seem to 
sharpen; see figures [7| and |H1 

If we wish to study the mobility of the solitons, moving breathers, etc the 




Figure 7. The energy (mass) of the breather for two coupled chains, as a function of 
the coupling constant Jo, and j3 = 2nm . From bottom to top, the curves correspond 
to: black g = 0, red g = 0.5, green g = 1, blue g = 1.2 and violet g = 1.5. In each case 
we have performed our calculations for both the increasing and decreasing J . 
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Figure 8. The profile of the DNLS breathers as a function of the coupling constant 
Jo, and (3 = 2nm~ 1 . The left figure shows g = 0, the right one g = 1.5. The sharpest 
curves occur for the smallest value of Jo, with a sequential broadening as Jo increases. 

standard technique j2H] in the field suggests the use of the smooth-profile (continuum) 
approximation, where we have added a term of the form 

- /A,i-i)/2 (15) 

to the left hand-side of the eigenvalue problem. Results of such calculations are given 
in figures 03 and EH Here the phase transition probably plays an important role as 
for small Jo we have a very localised solution which strongly violates the continuum 
approximation. In this case we are unable to find a moving soliton. For larger values of 
Jo we do find a smooth profile solution, which suggests that we have found the moving 
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Figure 9. The shape of the moving breather for two coupled chains, as a function of 
the coupling constant Jo, g — and (3 = 2nm _1 . The sharpest curves occur for the 
smallest value of Jo- 

soliton. This we expect to have an important impact on the charge transport properties 
of the system. Generally, the effect of the chirality of the molecule seems to sharpen the 
profile thus reducing the mobility! 
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Figure 10. The shape of the moving breather for two coupled chains, as a function 
of the coupling constant Jo, g = 1.5 and (3 = 2nm _1 . The sharpest curves occur for 
the smallest value of Jo broadening as Jo increases. 
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3.5. The molecule 

The model for the molecule can be specified in its most general form by solving a number 
of constraints, which we shall try to parametrise in as simple a way as possible. Looking 
at the figure, it seems easy to work with the following representation: (we use to 
denote the position where the ith base pair is attached to one of the chains, Sj for the 
other one) 

n = Zie z + ^dUi, Si = Zie z - ^du { , (16) 

with u a 2D unit vector in the xy plane. Hence it is natural to relate the angles 
9i introduced in (3) to the orientation of Ui in the xy plane. Our representation 
makes it possible to allow for stretching of the H bonds by adding a potential for d 
(as in the Peyrard-Bishop model). We shall use a molecular potential to implement 
an approximate constraint on the spacing in the strands of DNA, but it is not so 
straightforward to study the effect of bending of the double helix in this model, which 
could be quite important [2%] . 

The kinetic energy is relatively easy to evaluate 

\m r lt + s lt = \ M Yl ( Zie * + \ dv *) + { z ^ z - \ dUi 

i i ^ / t \ 

= M Y( Z l + \ d2u l) ■ ( 17 ) 

We can easily generalise this to a dynamical d; in this case the length of each rung 
becomes di(t), and the additional terms in the kinetic energy would be of the form \df t . 
However, we expect the effects of the variation of d to be small, so in this work we will 
keep d fixed. 

We shall now describe a set of simple potentials, all corresponding to properties of 
the bonds or the angles between bonds, that can be used to describe equilibrium DNA. 

If we insist on a molecule where the equilibrium shape is bent, such as DNA, we 
must impose an equilibrium angle between the bonds. One such potential would be 

U B oc Y [i( r i ~ x ~ r i)) 2 ~ q2 s] 2 

i 

+ x (s i+1 - Si)) 2 - a 2 s ] 2 . (18) 

Actually, this is too simple. In turning the outer product into a scalar, we have now 
allowed mixing of the left- and right-handed turns in the same molecule. Since DNA is a 
chiral molecule, we would like to add a term that prefers left-handed over right-handed 
twists. This is most easy to construct the continuum limit, where this translates into the 
requirement that the vector tangential to one of the strands points both upwards and 
into the plane, or in the opposite direction (downwards and out of the plane). Suppose 
that the point on the strand is specified by a vector d from the back bone (i.e., 8 lies 
in the xy plane), we write 8 = Se (in the notation used above, e = ±u, depending on 
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which strand we are on), with 

e = (cos sin 0,0)- (19) 

We can then write for the tangent vector 

t = a(- sin 9, cos 9,0) + 6(0,0,1), (20) 

and, to introduce chirality, we wish to add a bias towards a particular value of ab in our 
model. We can extract a by calculating e z • (e x t), and b — e z ■ t. Replacing t by the 
discrete analogue r i+1 — r*j, we find that 

+ u i+1 ) x ( (z i+1 - z,i)e z - ^{u i+1 - u, 



= ^e z ■ (uj+i x m) , (21) 
b r « e z • ^(zj+i - ^)e 2 - ^(ttj+i - t*i) 

= ^i+l — • (22) 

Here, due to the symmetry, a s = a r and b s = b r . A potential that achieves our aim is of 
the form 

U B = ^Qiiiizi+i - Zi)e z ■ (u i+ i x m)) , (23) 

i 

with g a function that has a minimum at /x, the chosen average parameter, and increases 
quickly as we move away from equilibrium. An example of such a function which in fact 
we have used in our simulations, is 

g(x) = b(x-fi) 2 . (24) 

We will also assume that each of the back-bones (the strands) can stretch; this will 
be built into our model by introducing a stretching potential of the bonds of the form 
(note that stretching for the other strand is constrained by the symmetry imposed by 
the straightness of our molecule) 

U s = J>(h-r m | 2 ) 

i 

d 2 

= ^2fi 2 ((zi+i - Zif + y(l - • «<)) , (25) 

i 

where / is has a minimum at the average distance / 2 , and increases quickly away from the 
minimum. As / describes a stretching potential we can either use a harmonic potential, 
f(x) = a(x — I 2 ) 2 , or a Morse potential, f(x) = a(e _ ^ _i ^ 1/2 — l) 2 . In our simulations 
we have used the simple quadratic form 

f{x) = a(x - I 2 ) 2 (26) 
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3. 6. A small simplification 

Before we proceed further we note that we can simplify the problem slightly by using a 
ID angle representation for u, 

u,i = (cos 6>j, sin^j, 0) . (27) 

We then find 

«M = # , (28) 

and 

U B = J29^(zi + i-z t )sm(6 t+1 -e i )) , (29) 

i 

as well as 

U s = X>((zi+i - z t ) 2 + y(l - cos(0 i+1 - 9A)) . (30) 

i 

Next we determine the equations of motion for all the fields of our system. They 
are given by: 

- hiijj itt = -2hwipi + ^2 Ji-j^j + \x*\%Vi + K<pi , (31) 
d 2 ■■ 

—M9i = + g'((z i+1 - Zi) sin(0 i+ i - 0A)( Zi+1 - zA cos(6» m - 9A 

- g'((zi - Zi-i) sm(9i - 0i-i))(zi - Zi-i) cos(6»i - 6>i_i) 

+ f'((z i+1 - zA 2 + d - (1 - cos(0 i+1 - sin(0 m - OA) 

- f'((zi-i - z t f + y (1 - cos(^ - i _ 1 )))^sin(^ - ^_0) - a9i 
+ terms depending on J (32) 

and 

2Mzj = + g'((z i+1 - ^) sin(fl m - 9i)) sin(0 m - 0*) 

- g'((zi - sin(6»i - 6>;_i)) sin(0j - 6>i_i) 

+ f((z m - ^) 2 + — (1 - cos(0 m - 0i)))2(^ +1 - ^) 
d 2 

- /'((zi - Zi_i) 2 + — (1 - cos(6»i - 6>j_i)))2(^ - Zi_i) - aij 

+ terms depending on J . (33) 

Please note that we have added a damping term to both the z and 9 equations (the term 
proportional to a). This helps us to lower the energy in our system and to determine 
its lowest energy states which will be solitonic if the system possesses such solutions. 

The only remaining question is the choice of g. A sensible choice seems to be: 
9n( x ) — K x ~ I 1 ) 2 an d this is the expression we used in our simulations. 
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3. 7. J dependence 

Next we add the J dependent terms which will control electron's mobility. Although, 
in principle, J can depend on many variables two obvious dependences first spring to 
mind: 

(i) The tunnelling rate is a function of the shortest distance between two spatial points, 
i.e., r {j = \r { - Vj\, or 

(ii) The tunnelling rate depends on the distance along the DNA backbone (the sucrose 
molecules), fy = YjkZj \ r i ~ r i+i\- 

The first case is obviously easier to implement, but shall not be used in this work, 

The second case is more complicated. Nonetheless; we can write Jy = h(Rij) and 
find that the additional terms in the EoM for 9, Eq. ()32j) are of the form 

d 2 sin(0 i+ i - 6i) , 
> h (-Rjw)Vw + {k^l) 

k<i,l>i 

_ d^n(9i-i-0i) J2 h'^^ + ^^i) (34) 

r i-i i 

1 i '' k<i,l>i 



and for z (Eq. 

d 2 Zi(Zi — 



r * +1 ' 4 k<i,l>i 



h'(R kl )r k ^i + (k^l) 



f/J:,(: " E h'(R kl )m+(k~i) ■ (35) 



r i~i % 

k<i,l>i 



3.8. Choice of parameters 

Next we attempt to determine the realistic (i.e. of some relevance in biophysics) range 
of our parameters. First we look at the static solutions of the "molecular" part of the 
energy, and minimise both the bending and the stretching potentials; we find that 



(z i+ i - Zi)e z ■ (ui x Uj+i) = /i , (36) 

d? d 2 

(z i+t - Zi) 2 - — u i+l ^ = I 2 - — . (37) 

Assuming uniformity we get 

sin 56* = /i , (38) 

5z 2 + d 2 sin 2 (59/2) = l 2 , (39) 

where 5z = (z i+ i — z^), 56 is the angle between two subsequent it's. Thus it is simple to 

find a value for /i and / given 5z, d and 56, which we take from the studies of the real 
DNA. The standard values are j2Zj 

d = 2.39 nm , (40) 

5z = 3.38/10 = 0.338 nm , (41) 

56 = 2tt/10 = 0.628 . (42) 
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This gives us 



0.20 nm 



I = 0.81 nm 



(43) 



If we look at all the minima for these values of /z, I and d we find that we have a second 
solution 



this indicates an obvious weakness of our model — it will lead to coexisting shapes! 
4. Numerical Studies 

We have performed several numerical studies of our system. We have taken sensible 
parameters; i.e. those that appeared not to be completely inconsistent with the natural 
physical constraints but, at the same time, are quite useful to gain some understanding 
of the nature of the solutions. We have taken h — M — 1, u — 1 a = 6, Z 2 = 0.7 
(J2HJ), b = 6,/i = 0.2 (J21j), K — 0.1 (jSJ). The damping for the molecular motion was set 
at a = 0.02. We also varied the value of J from 0.01 to 10. 

We have found that for small values of J the system behaved in a sensible way - 
i.e., the molecules were quite rigid and the electron field modified their behaviour very 
little. The electron field was solitonic in nature and it affected the molecule primarily 
in the places where it was nonzero. For larger values of J — i.e., J > 0.035 the electron 
field dispersed and its effects were more pronounced. This can be most easily seen by 
looking at the figures which we include. To obtain them we used a 101 "rung" model of 
the twisted DNA-ladder, i.e., with 10 complete turns. We note that the transition point 
from the localised to the delocalised regime appears to have moved a bit as compared 
to the DNLS, see below. 

In our first simulations we looked at J = 0.02. As can be seen from figure ^2 m 
this case, the kinetic energy decayed rapidly, and so did the potential, leaving us with 
a breather soliton, see figure H21 

In this figure we show that the amplitude (square-root of the probability) of the 
electron field at one rung (i.e. summed over the two sites connected by that rung) 
performs semi-regular oscillations. The simulation was started out with all probability 
on one site of one backbone. We note also that some probability starts propagating 
outwards. This is then reflected from the boundaries, returns to the soliton and makes 
the breather's oscillations somewhat less periodic resulting in a small sea of background 
oscillations. 

Some aspects of the periodicity can also be seen in figure ITB1 where we have plotted 
the time dependence of the total probability on one of the two backbones (i.e. summed 
over all 101 sites on that backbone). We observe a fast decay to 0.5, and then small 
oscillations around this value, related to the periodicity of the breather. 

We note that the presence of the soliton leads to a pronounced deviation of 9 and z 
from their equilibrium values, with a sharp jump at the position of the soliton ffigurelT4*j). 
There seems to be a stronger correlation between these two quantities than one would 
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Figure 11. Energetics of the relaxation of a 101-rung chain for J = 0.02. The upper 
panel with inset shows the decrease of the mechanical energy, and the lower panel the 
same for the total and electronic energy 



have naively expected, but the overall picture shows that there is a deformation and 
relaxation of the double helix. 

Next we have studied a transitional case, namely, J = 0.03. 

Once again the relaxation to the equilibrium has been quick and uneventful, at 
least as far as the energies are concerned, see figure ITB1 

The observed time evolution of the amplitude shows that we do have a stationary 
solution, but the effect of the small part of the probability density that propagates 
is much more pronounced, and makes things look much more chaotic (figure IT6|) . In 
the early time picture we can also see that these probability waves get reflected back 
towards the soliton (and also get reflected by the soliton). Clearly, the interesting 
dynamics observed here deserves further study. 

We also note similar structure in the oscillations from one backbone to the next, 
see figure El The oscillations are more pronounced, and their envelope is probably due 
to the probability bouncing from the boundaries. The fastest oscillation is then, almost 
certainly, the frequency of the breather. 

The sharp jump of z and 9 at the soliton position can again be seen in figure 

The situation changes drastically if we increase J to 0.04. We still have a well 
defined relaxation, figure E3 but there is no longer a stable soliton. 
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Figure 12. Time evolution of the occupation as a function of the position along the 
backbone (sum over both) for J = 0.02. t gives a time slice number (100 slices used 
for each figure). The top figure shows the early time evolution, the lower one the late 
time steady state breather solution. A dynamical representation of the evolution is 
available at http://walet.phy.umist.ac.uk/DNA/Jp02.php 



In figure EH we demonstrate this. We note a fast decay of the soliton, even before 
the reflected waves have a chance to return to it. This suggest the absence of a stable 
breather. At later times, as shown in the lower panel, we see that a soliton-like structure 
gets reassembled at a different position. This may well indicate that at this point the 
solitons are only marginally unstable. 

The envelope of the norm also shows more structure, even though this structure 
seems to be overlaid on a very regular pattern, figure HT\ 

The coupling to the position of the mechanical structure of the double helix is 
substantial, as can be seen from the results shown in figure E21 
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Figure 13. Time evolution of the total occupation of one of the backbones as a 
function of time for J = 0.02. The insets show both early and late oscillations of this 
quantity. 
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Figure 14. Deviation of 8 (red squares) and z (black circles) from their equilibrium 
values for the final simulation point, t = 19318, as a function of the lattice point 
labelled by i. The coupling parameter J = 0.02. 



5. Conclusions 

We have constructed a model that describes the chiral movement of electrons in a double 
helix. We have shown how it describes both localised and delocalised charge, depending 
on parameters. 

We have shown that the dynamics of a chiral molecule chain can be quite involved 
and interesting, leading to a complicated interplay of stretching and twisting modes 
of the underlying double helix generated by the movement of the electrons. As Jo 
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Figure 15. Energetics of the relaxation of a 101-rung chain for J = 0.03. The 
upper panel with inset shows the decrease of the mechanical energy, the lower panel 
the decrease of the total and the electronic energy. 



increases the system exhibits a sharp transition from localised solitons (charge density) 
to a completely delocalised system with a complicated temporal structure. In each case 
we see an anticorrelation between z and 9, with the soliton deforming the lattice by 
locally reducing the spacing of the double helix and increasing its twist. 

Our next step will be to try to extend these calculations to charge transport, where 
electrons come in on one backbone at one end, and leave through the same or the other 
backbone at the other side. 
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Figure 21. Time evolution of the total occupation of one of the backbones as a 
function of time for J = 0.04. 




Figure 22. Deviation of 9 and z from their equilibrium values for the final simulation 
time, t = 13339.5, as a function of the lattice point label i and J = 0.04. 



